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Abstract. We investigate and contrast, via entropic sampling based on the Wang- 
Landau algorithm, the effects of quenched bond randomness on the critical behavior of 
two Ising spin models in 2D. The random bond version of the superantiferromagnetic 
(SAF) square model with nearest- and next-nearest-neighbor competing interactions 
and the corresponding version of the simple Ising model are studied and their general 
universality aspects are inspected by a detailed finite-size scaling (FSS) analysis. We 
find that, the random bond SAF model obeys weak universality, hyperscaling, and 
exhibits a strong saturating behavior of the specific heat due to the competing nature 
of interactions. On the other hand, for the random Ising model we encounter some 
difficulties for a definite discrimination between the two well-known scenarios of the 
logarithmic corrections versus the weak universality. Yet, a careful FSS analysis of our 
data favors the field-theoretically predicted logarithmic corrections. 
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1. Introduction 

The understanding of the role played by impurities on the nature of phase transitions 
is of great importance, both from experimental and theoretical perspectives. First- 
order phase transitions are known to be dramatically softened under the presence of 
quenched randomness [H El El HI El E], while continuous transitions may have their 
exponents altered under random fields or random bonds [21 [7J E]. There are some 
very useful phenomenological arguments and some, perturbative in nature, theoretical 
results, pertaining to the occurrence and nature of phase transitions under the presence 
of quenched randomness 0, El QUI EEH DEI E3j. The most celebrated criterion is that 
suggested by Harris [7J. This criterion relates directly the persistence, under random 
bonds, of the non random behavior to the specific heat exponent a p of the pure system. 
According to this criterion, if a p is positive, then the disorder will be relevant, i.e., under 
the effect of the disorder, the system will reach a new critical behavior. Otherwise, if a p 
is negative, disorder is irrelevant and the critical behavior will not change. The value 
a p = is an inconclusive, marginal case. The 2D Ising model falls into this category, 
it is the most studied case, but is still controversial |14j . In general and despite the 
intense efforts of the last years on several different models, our current understanding 
of the quenched randomness effects is rather limited and the situation appears still 
unclear for both cases of first- and second-order phase transitions. The present paper 
aims to contribute to our knowledge of the effects of quenched bond disorder on second- 
order phase transitions. In this quite active field of research, the resort to large scale 
Monte Carlo simulations is often necessary and useful, and following this recipe, we 
applied some recently developed by our group numerical schemes - see reference [15] 
and references therein - on two types of random bond Ising models. 

In the first part of the paper we present an extension of our numerical investigation 
of a random bond spin model in 2D with competing interactions, known as the random 
bond square SAF model [16] . Details of the pure and random version of this model with 
competing interactions and the motivation of such a study will be given below in the 
corresponding Section, but can also be found in our recently published Letter [16J and 
in other related papers in the literature [17j. Furthermore, in the second part of our 
work, by a parallel study - using the same numerical techniques - we attempt to shed 
new light into the well-known random bond version of the 2D (simple) Ising model. 
Our investigation will be related to the extensive relevant literature concerning this 
case [HI HH [201 1211 [221 [23l [2H ESI [26l [271 (211 [291 EQl EH [32l [33l EH ESI [361 EH EH EH1 

sniaiiaiaaasiaaiia in particular, our 

discussion will focus on the main point of the last two decades, concerning the two well- 
known conflicting scenarios, namely the logarithmic corrections [181 EH EU E2] versus 
the weak universality scenario [HTJ E21 EHl [60] . 

The rest of the paper is organized as follows: In the next Section we outline an 
extensive entropic sampling program. This program is based on (i) the Wang-Landau 
(WL) method [61], (ii) the dominant energy restriction scheme [62J, and (iii) a second 
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stage improvement that combines the WL method [HI] and some new ideas [T51 [T6"| 63j. 
suitable for the study also of disordered systems. Our FSS analysis of the numerical 
data and the corresponding discussion of the random bond versions of the square SAF 
model and the 2D Ising model are presented in Sections [3] and HI respectively. Finally, 
our conclusions are summarized in Section [5j 

2. Simulation schemes and numerical details 

Importance sampling methods have been for many years the main tools in condensed 
matter physics and critical phenomena [64j [65l [66l [67J [68] • However, for complex 
systems, effective potentials may have a rugged landscape, that becomes more 
pronounced with increasing system size. In such cases, these traditional methods 
become inefficient, since they cannot overcome such barriers in the state space. A 
large number of generalized ensemble methods have been proposed to overcome such 

problems [SH EZl IEHI IBSl ISIl EH ISl ISl ISl ISl ISl 1^1 1^1 IHOl ISH IH2|- One 
important class of these methods emphasizes the idea of directly sampling the energy 
density of states (DOS) and may be called entropic sampling methods [67]. In entropic 
sampling, instead of sampling microstates with probability proportional to e~ l3E , we 
sample microstates with probability proportional to [G(E)]^ 1 , where G(E) is the DOS, 
thus producing a flat energy histogram. The prerequisite for the implementation 
of the method, is the DOS information of the system, a problem that can now be 
handled in many adequate ways via a number of interesting approaches proposed in 
the last two decades. The most remarkable examples are the Lee entropic [69, [70], the 
multicanonical [73~1I7I]. the broad histogram [71], the transition matrix [72], the WL [61J, 
and the optimal ensemble methods [82] . 

In particular, the WL algorithm [61] is one of the most refreshing variations of 
the Monte Carlo simulation methods introduced in the last years. The algorithm has 
already been successfully used in many problems of statistical physics, biophysics, and 
others [831 Ell ESI EHl [871 EH1 EH |90l EH E21 [931 EH [951 [961 EZl EH]- To apply the WL 
algorithm, an appropriate energy range of interest has to be identified and a WL random 
walk is performed in this energy subspace. Trials from a spin state with energy E^ to a 
spin state with energy Ef, using local spin flip dynamics, are accepted according to the 
transition probability 

p(Ei — > Ef) = min 

During the WL process the DOS G(E) is modified (G(E) -> / * G(E)) after each 
spin flip trial by a modification factor / > 1. In the WL process (j = 1,2,..., jj) 
successive refinements of the DOS are achieved by decreasing the modification factor 
fj. Most implementations use an initial modification factor fj = i = e ~ 2.71828. . ., a 
rule fj + i = J~fji an d a 5% — 10% flatness criterion (on the energy histogram) in order 
to move to the next refinement level (j — > j + 1) [61]. The process is terminated in a 
sufficiently high-level (/ ~ 1, whereas the detailed balanced condition limit is / — > 1). 



G{Ej) 
G(E f ) 



(1) 
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In the last few years, there have been several papers dealing with improvements and 
sophisticated implementations of the WL iterative process [63j EH EH EE EH EH HOOj 
10 lj. The present authors have introduced a dominant energy subspace implementation 
of the above entropic methods, called critical minimum energy subspace (CrMES) 
method [62] . This is a method of a systematic restriction of the energy space, with 
increasing lattice size, by which one can determine all finite-size thermal anomalies of the 
system from the final accurate DOS, and also other (magnetic) anomalies of the system 
by accumulating appropriate histogram data in the final almost entropic stage of the 
process. The (WL) random walk takes place in a restricted energy subspace (E%, E 2 ) and 
this practice produces an immense speed up, without introducing observable errors. It 
has been shown that the method can be followed successfully in pure systems undergoing 
second- or first-order phase transitions [62l [85] and also in more complex systems with 
complicated free energy landscapes, such as the 3D random field Ising model [15] . 

For the simulation of the random bond models considered in this paper, we followed 
the general framework of the implementation of the above described scheme. In 
particular, we followed most of the details of the implementation applied recently to 
the 3D random field Ising model and outlined for the random bond version of the SAF 
model in our recently published Letter [16]. In these papers, a two stage strategy has 
been followed. In the first stage, a multi-range (multi-R) WL method was applied, where 
the total energy range was split in many subintervals [61] and the DOS's of these separate 
pieces were then joined at the end of the process. The WL refinement levels used in 
this first multi-R WL stage (j = 1, . . . , j«), were as follows: ji = 18 for L < 80, ji = 19 
for 80 < L < 120, and ji = 20 for L > 120. In the second stage of the simulation 
(WL refinement levels: j = ji — j/), a more demanding multi-R - but with larger 
energy pieces - or an one-range (one-R) approach was carried out. The identification 
of the appropriate energy subspace (Ei,E 2 ) for the entropic sampling of each disorder 
realization was carried out by applying our CrMES restriction [62] and taking the union 
subspace at both pseudocritical temperatures of the specific heat and susceptibility. This 
union subspace, extended by 10% from each side, low- and high-energy side, is in most 
cases sufficient for an accurate estimation of all finite-size anomalies. The identification 
of the appropriate energy subspace was carried out in the first multi-R WL stage, using 
originally a very wide energy subspace. After the first identification, the same first stage 
process (j = 1, . . . , ji) was repeated several times, typically ~ 4 — 6 times, in the new 
restricted energy subspace. From our experience, this repeated application of this first 
stage multi-R WL approach greatly improves accuracy, and then the resulting accurate 
DOS is used again for a final and more accurate redefinition of the subspace (Ei,E2), 
in which the final entropic scheme (second stage) is applied. Thus, the final entropic 
WL stage (j = j, — j/), where jf = ji + 4 in all cases, was carried out in this accurately 
defined subspace in an one-R or in a multi-R approach. For the present models, it 
was found that a final multi-R approach with large subranges (see below) is in fact 
sufficiently accurate. Therefore, since the multi-R of the original WL scheme improves 
efficiency, we applied, for most of our simulations, this multi-R WL approach also in 



Quenched bond randomness in marginal and non-marginal Ising spin models in 2D 5 

the final entropic stage, using three times larger energy subintervals than in the initial 
multi-R stage. The energy subintervals of the first stage where chosen to correspond to 
rather large subspaces, with their sizes depending on the disorder strength. Taking the 
pure system as reference, these energy subintervals could be chosen of the order of 50 
to 100 energy levels, depending on the lattice size. In the disorder case, the subinterval 
sizes are multiplied by the factor induced due to the new multiplicity of energy levels, 
giving for instance, a factor 4 for the disorder strength r = 3/5 = 0.6, where r is the 
ratio of weak over strong bond interactions (see also next Section). 

As pointed out above, the need of using in the final stage the described multi-R 
approach, instead of an one-R approach, is a consequence of the slow convergence at 
the high WL levels. It is possible to overcome this slow convergence by using a looser 
flatness criterion or an alternative Lee entropic final stage, as proposed in reference [70] 
and applied by the present authors [THJ - However, recently a different alternative has 
been proposed by Belardinelli and Pereyra (BP) [63], which is free of the application 
of the energy-histogram flatness criterion. Following their proposal, one is using, in the 
final stage, an almost continuously changing modification factor adjusted according to 
the rule In/ ~ t~ x . Since t is the Monte Carlo time, using a time-step conveniently 
defined proportional to the size of the energy subinterval, the efficiency of this scheme 
is independent of the size of the subintervals and therefore the method provides the 
same efficiency in both multi-R and one-R approaches. Furthermore, from the tests 
performed by these authors, and also from our comparative studies in the 2D pure 
Ising model (unpublished), the error-behavior of this method seems superior to the 
original WL process, improving to some extent the saturation-error problem of the 
WL method. Accordingly, we have also applied this alternative route for the final 
stage of our simulations using an one-R approach. In particular, the disorder strength 
case r = 9/11 = 0.818 of the random bond square SAF model, and all simulations 
corresponding to the larger sizes L = 160 and L = 200 of the random bond Ising model 
were carried out using this alternative. From our comparative tests, we found that both 
approaches (the BP and the multi-R WL approach) produced very accurate results, 
with the BP approach giving superior estimates for the pseudocritical temperatures of 
the models. 

Both disordered models were simulated for two values of the disorder strength r. 
For the random bond square SAF model we chose the values r = 9/11 = 0.818 and 
r = 3/5 = 0.6, whereas for the random bond Ising model the values r = 3/5 = 0.6 
and r = 1/7 = 0.142 have been considered. Square lattices, using periodic boundary 
conditions, with linear sizes L in the range L = 20 — 120 or L = 20 — 200 (disorder 
strength case r = 1/7 = 0.142) were used. A number of 200, disorder realizations 
was generated and simulated for each disorder case and lattice size. Even for the larger 
lattice sizes the statistical errors of the WL method (WL-errors), used for the estimation 
of thermal and magnetic properties of a particular realization, were found to be much 
smaller than the statistical errors of the disorder averaging, coming from the fact that we 
have used a finite number of 200 disorder realizations. Therefore, the WL-errors are not 
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shown in our graphs, whereas the latter errors of finite disorder sampling are presented 
in our figures as error bars. The mean values over disorder will be denoted as [. . .] av , the 
corresponding maxima as [. . .]*„, and finally the individual maxima as [. . .*} av - Since 
in the fitting attempts of the following Sections, we have used mainly data from the 
peaks of the disorder averaged curves (i.e. [C]* v ), their finite disorder sampling errors 
are the relevant statistical errors to be used in the fitting attempts. These errors have 
been estimated by two similar methods, using groups of 25 to 50 realizations for each 
lattice size and the jackknife method or a straightforward variance calculation (blocking 
method) [67]. The jackknife method yielded some reasonably conservative errors, about 
10 — 20% larger than the corresponding calculated standard deviations, and are shown 
as error bars in our figures. Finally, let us point out that in all cases studied, the sample- 
to-sample fluctuations for the individual maxima are much large than the corresponding 
finite disorder sampling errors. 

3. The random bond square SAF model: Competing interactions 

In this Section we extend our investigation on the effects of quenched bond randomness 
on the square Ising model with nearest- (J nn ) and next-nearest-neighbor (J nnn ) 
antiferromagnetic interactions. In zero field, the pure square SAF model, is governed 
by the Hamiltonian: 

7~tp Jnn ^} ] SiSj + Jnn^SiSj, (2) 

<iJ> 

where here both nearest- (J nn ) and next-nearest-neighbor (J nnn ) interactions are 
assumed to be positive. It is well-known that the model develops at low temperatures 
SAF order for R = J nn /J nnn > 0.5 [1021 11031 1104] and by symmetry the critical 
behavior associated with the SAF ordering is the same under J nn — ► —J nn . For 
the case R — 1, that we deal with, the pure system undergoes a second-order 
phase transition, in accordance with the commonly accepted scenario for many years 
of a non-universal critical behavior with exponents depending on the coupling ratio 
R [1021 HD31 [1051 M [107]. The recent numerical study of Malakis et al [108] has 
refined earlier estimates [1031 1106] for the correlation length exponent v and values 
very close to those of the 2D three-state Potts model z/ p (Potts)= 5/6 [109] were 
obtained. From the FSS of the pseudocritical temperatures [108] it was found that 
u p (SAF;R = 1)= 0.8330(30) and the subsequent study of Monroe and Kim [110j . 
using the Fisher zeroes of the partition function, yielded a quite matching estimate: 
z/ p (SAF;i? = 1)= 0.848(1). Furthermore, from the FSS of the specific heat data an 
estimate for the ratio ct p /v p = 0.412(5) was also found [108] . Finally, from the magnetic 
data and in accordance with an earlier conjecture of Binder and Landau [103] . Malakis 
et al [108J found additional evidence of the weak universality scenario [591 EQ] and 
obtained the values (3 p /v p = 0.125 and 7 P /^ P = 1.75. The values of the above three 
ratios of exponents satisfy the Rushbrook relation, assuming that v p = 0.8292, which is 
very close to the estimate obtained from the shift behavior of the SAF R = 1 model, thus 
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Figure 1. Size dependence of the maxima of the specific heat for the pure (filled 
squares; data taken from reference [108] ) and the random bond (filled and open circles 
for r — 0.818 and open triangles for r = 0.6) square SAF model. The inset shows a 
power law fit for the case r = 0.818 for L > 80 giving a negative value for the exponent 
ajv of the order of —0.12(6). 



providing self-consistency to the estimation scheme. From these results, it is tempting 
to conjecture, as was pointed out in reference P3], that the SAF model with R = 1 obeys 
the same thermal exponents with the 2D three-state Potts model [y v = 5/6 = 0.833 . . . 
and a p = 1/3 = 0.333... |109j ). but the respective values of the magnetic critical 
exponents are different [fl P lv v = 2/15 = 0.133 . . . and j p /u p = 26/15 = 1.733 ... for the 
2D three-state Potts model [T09] l 

Considering now the random bond distribution [16] 

p(J ij ) = \[5(J ij -J 1 ) + 5(j ij -j 2 )\; ^r^ = l; r = ^, (3) 

where the ratio r stands for the disorder strength, the Hamiltonian of equation (TjQ) is 
transformed into the following disordered one 

ft = JijSiSj + J2 JijSiSj- ( 4 ) 

The critical behavior of the above defined random model for the case r = 3/5 = 0.6 has 
been outlined in reference |16j . where apart from the verification of the weak universality 
scenario, a strong saturating behavior of the specific heat has been witnessed. Here, we 
extend our study to a weaker disorder strength, namely the case r = 9/11 = 0.818. 

Let us start the presentation of our results with the most striking effect of the 
bond randomness on the specific heat of the square SAF model. In figure [1] we contrast 
the size dependence of the specific heat maxima of the pure (filled squares) and the 
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Figure 2. (a) Simultaneous fittings of several pseudocritical temperatures of the 
random bond square SAF model for the two values of r considered: r — 0.818 and 
r = 0.6. (b) Log-log plot of the maxima of the average second-order logarithmic 
derivative of the order-parameter of the random bond square SAF model also for 
r = 0.818 and r — 0.6. The inset shows a log-log plot of the maxima of the average 
logarithmic derivatives of the order-parameter of first-, second-, and fourth-order, for 
the case r = 0.818. Linear fits are applied for L > 30. 

random bond model (filled and open circles for r = 0.818 and open triangles for the 
case r = 0.6). For the case r = 0.818 we show two data set points for the specific 
heat, corresponding to the two averaging processes discussed previously in Section EJ 
Note that the error bars for the quantity [C*] a „ shown reflect the sample-to-sample 
fluctuations, whereas all other error bars are statistical errors due to the finite number 
of disorder realizations (jackknife errors discussed in Section [2]). The suppression of the 
specific heat maxima is clear for both disorder strength values and of course it is much 
stronger for the case r = 0.6, for which a clear saturation is observed even for the smaller 
sizes shown (L = 40). For the present value r = 0.818 we show in the inset of figure [T]a 
power law fitting attempt of the form [C]* av ~ + bL a l v for sizes L > 80 which gives 
a negative value for the exponent ajv of the order of —0.12(6). Notably, this value of 
ajv will be shown to be compatible with the one obtained by an alternative method via 
the Rushbrook relation. 

In figure [2](a) we present the FSS behavior of several pseudocritical temperatures of 
the model T\z\* , i.e. the temperatures corresponding to the maxima of several average 
curves, such as the specific heat {Z = C), the magnetic susceptibility [Z = x), the 
absolute order-parameter derivative with respect to the temperature [Z = 9< ^ > ), and 
the logarithmic derivatives of several powers of the order-parameter with respect to the 
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Figure 3. (a) Log-log plot of the size dependence of the average magnetization of the 
random bond square SAF model at the estimated corresponding critical temperatures 
for r — 0.818 and r = 0.6. (b) Log-lop plot of the size dependence of the maxima of the 
average magnetic susceptibility of the random bond square SAF model for r = 0.818 
and r = 0.6. In both panels, the solid and dotted lines are corresponding linear fits for 
L > 30. 

temperature [Z = din<M n > ^ n = and 4). For the case r = 0.818 all the above 
6 mentioned pseudocritical temperatures are presented and are compared to the case 
r = 0.6 [16J. The lines show two sets of simultaneous fitting attempts, according to the 
shift relation 

T [z] , v =T c + bL-V\ (5) 

giving T c = 2.058(8) for r = 0.818 and T c = 1.980(9) for r = 0.6, respectively, for 
the critical temperature of the disordered model. Note that the corresponding critical 
temperature of the pure system is T c . p = 2.0823(17) [108] . The values for x 2 /DoF of 
the fits shown depend on the method used to evaluate the statistical errors (jackknife 
or simple standard deviation errors) and in all cases studied in this paper vary in the 
range 0.1 — 0.5. A first estimation of the critical exponent v of the correlation length 
is obtained from the above shift behavior and is v — 1.050(13) and v = 1.080(20) for 
r = 0.818 and r = 0.6 respectively, as illustrated in the graph. An alternative estimation 
of the exponent v is attempted now from the FSS analysis of the logarithmic derivatives 
of the order-parameter [111] with respect to the temperature [H 1112] 
01n(M") (M n E) 

which scale as L x l v with the system size. In figure |2]^b) we consider in a log-log scale 
the size dependence of the maxima of the average second-order logarithmic derivative 
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of the order-parameter for r = 0.818 (open circles) and r = 0.6 (open triangles). The 
solid and dotted lines are corresponding linear fits whose slopes provides respectively 
estimates for 1/v and thus for the exponent v. It is clear for the figure that the slopes of 
the lines are different and the results we obtain from the linear fits are are v = 1.047(10) 
and v = 1.090(12) for r = 0.818 and r = 0.6, respectively. In the corresponding inset of 
figure |2fb) we present the first- (filled squares), second- (open circles), and fourth-order 
(open triangles) maxima of the average over the ensemble of realizations logarithmic 
derivatives of the order-parameter for the case r = 0.818. The solid, dashed, and dotted 
lines shown are corresponding linear fits whose slopes provide an average estimate for 
the exponent v of the order of v = 1.046(10). These results for the exponent v for both 
values of the disorder strength r compare favorably with the estimation of v from the 
above shift behavior shown in panel (a) of figure O Thus, in comparison with its value of 
the pure model, the exponent v for the disordered model shows an increase of the order 
of 20% for the case r = 0.818 and 30% for the case r = 0.6, reflecting, in both cases, the 
strong influence of the disorder on the thermal properties of the system. Noteworthy 
that, our estimates for the exponent v are in agreement with the inequality v > 2/D 
derived by Chayes et al [8] for disordered systems. 



Turning now to the magnetic properties of the model we present in figure [3](a) in 
a log- log scale the FSS behavior of the average order-parameter for r = 0.818 (open 
circles) at T c = 2.058 and r = 0.6 (open triangles) at T c = 1.98. The straight 
lines show linear fits for L > 30 with a slope of 0.124(6) for r = 0.818 and 0.126(5) 
for r = 0.6. Thus, these two estimations for the ratio (3/v indicate that although 
the exponent f3 increases with disorder, the ratio (5/v remains unchanged to its pure 
value, i.e. f3/v = fi v jv v = 0.125 [TB I I108j . Correspondingly, we show in panel (b) of 
figure [3] the behavior of the magnetic susceptibility that provides estimates for the ratio 
7/V. The open circles refer to the case r = 0.818 and the open triangles to the case 
r = 0.6. The solid and dotted lines are linear fits giving the values 7/1/ = 1.749(9) 
and 7/// = 1.751(10), respectively. Thus, the ratio 7/1/ maintains the value of the 
pure model for both disorder strength values considered. The ratios (3/v and j/u for 
the disordered square SAF model appear to be the same with the corresponding ratios 
of the pure square SAF model but different from those of the 2D three-state Potts 
model. Therefore, our results reinforce both the weak universality scenario for the pure 
SAF model, as first predicted by Binder and Landau |103j . as well as the generalized 
statement of weak universality in the presence of bond randomness, given by Kim [113] 
and concerning also the 2D random bond three-state Potts ferromagnet. 

Finally, having estimated the ratios (3/u, 7/1/, and the exponent u, we attempt to 
estimate the specific heat exponent a using either the Rushbrook (a + 2(3 + 7 = 2) 
or equivalently, since 2(3 /v + 7/1/ = 2, the hyperscaling (2 — a = Dv) relation. For 
the case r = 0.818 we find a value ajv = —0.09(4) which is in agreement with 
the value ajv = —0.12(6) estimated from the fitting of the larger specific heat data 
shown in the inset of figure [TJ In the case of a non-divergent specific heat it is quite 
interesting to consider also, rather than the specific heat, the internal energy scaling at 
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the estimated critical temperature. This provides an estimate for the exponent ratio 
(a — 1)/V, which may be more precise [114] . We finally carried out such a fitting using 
the larger sizes (L > 80) on the values of the critical energy (at T c = 2.058) and we found 
(a — l)/u= —1.04(4). This result, when combined with the estimate for v from the shift 
behavior [y = 1.050(13)) gives a value a/v — —0.09(3), an intriguing coincidence, in 
full agreement with the above value obtained for this ratio, via the Rushbrook relation. 
For the case r = 0.6 [16] , the saturation effect is much more stronger and an even more 
negative value for the specific heat exponent a is obtained a = —0.17(4). From the above 
and in agreement with our preliminary observation in reference [16], it turns out that this 
strong saturating behavior of the specific heat is completely different from the behavior 
of the 2D random bond three-state Potts ferromagnet |113l 1115] . There, a specific heat 
diverging behavior is obtained for disorder strengths r = 0.9, 0.5, and 0.25 |113j and 
an increasing but progressively saturating behavior is obtained only for the very strong 
disorder r = 0.1 [115j . This behavior of the specific heat of the random bond square 
SAF model may presumably be attributed to the competitive nature of interactions, 
responsible for the observed sensitivity of the SAF model to bond randomness. 

4. The marginal case of the 2D random bond Ising model 

The Hamiltonian of the random bond version of the 2D Ising model system is given by 



<i,j> 

where again the implementation of the bond disorder follows the binary distribution (J3]) 
of ferromagnetic interaction strengths. With this distribution the random Ising system 
exhibits a unique advantage that its critical temperature T c is exactly known |116[ 1117] 
as a function of the disorder strength r through 



where r = J2/J1 and k B = 1. This gives the opportunity of carrying out the FSS 
analysis at the exact T c (r), a practice that highly reduces statistical errors. Furthermore, 
one may check for accuracy his numerical scheme by comparing the estimated critical 
temperature (via the shift behavior) with the exact result, as we have also done here. 

At this is point, it is useful to briefly discuss the two main and mutually excluded 
scenarios [H] mentioned in the introduction. The logarithmic corrections scenario is 
based on the quantum field theory results of Dotsenko and Dotsenko [IS] , and of 
Shalaev [19], Shankar [21], and Ludwig [22J. According to this scenario - supported 
theoretically for the case of weak disorder - the presence of quenched disorder changes 
the critical properties of the system only through a set of logarithmic corrections to 
the pure system behavior. The specific heat C is expected to diverge on approach 
to the critical temperature T c in a double logarithmic form: C oc ln(lnt), where 
t — \T — T c \/T c is the reduced critical temperature. On the other hand, according 
to the weak universality scenario [3T|, [32l [591 EQ] , critical quantities, such as the zero 



n 




(7) 



sinh (2 Jx/T c ) sinh {2rJ x /T c ) = 1, 



(8) 
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Figure 4. Critical behavior of the 2D random bond Ising model for r = 3/5 = 0.6. 
(a) Simultaneous fitting of two pseudocritical temperatures defined in the text and 
estimation of the correlation length exponent, (b) FSS behavior of the averaged specific 
heat at T c (r = 0.6). The solid and dotted lines are corresponding double logarithmic 
and power law fits (see equations © and (fl0|) in the text). The inset shows the specific 
heat data as a function of the double logarithm of the lattice size. A linear fit for L > 60 
is applied, (c) Log-log plot of the averaged magnetization at T c (r = 0.6). (d) Log-log 
plot of the averaged magnetic susceptibility at T c (r — 0.6). In both panels (c) and (d) 
a linear fit is applied. Note that all fitting attempts shown have been performed for 
L > 20. 

field susceptibility, magnetization, and correlation length display power law singularities, 
with the corresponding exponents 7, (3, and v changing continuously with the disorder 
strength; however this variation is such that the ratios 7/z/ and f3/v remain constant 
at the pure system's value. The specific heat of the disordered system is, in this case, 
expected to saturate [3"T] . 

Two significantly different values of the disorder strength, namely the cases r = 
3/5 = 0.6 and r = 1/7 = 0.142 have been investigated by our numerical scheme 
and our data will be presented and analyzed below. Due to marginality, the case 
r = 0.6, that produced in the SAF model a 5% temperature decline, gives here, 
through equation (JHI) a critical temperature T c (r = 0.6) = 2.22419 . . ., very close 
(2%) to the corresponding critical temperature of the pure system (T c . p = 2.26918. . .). 
Thus, we may call the first case (r = 3/5 = 0.6) a weak disorder and the second case 
(r = 1/7 = 0.142) a strong disorder case, actually producing a significant temperature 
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decline T c (r = 0.142) = 1.77910.... Our simulations are extended to lattice sizes in 
the range L = 20 — 120 for the weak disorder and in the range L = 20 — 200 for the 
strong disorder, in the hope that we will observe the true asymptotic behavior of the 
model. For the calculation of the statistical errors due to the finite number of simulated 
realizations, we have followed again the practice outlined in Section El The values of 
X 2 /DoF of all the fits shown below are again in the range 0.1 — 0.5. 

We start the presentation of our results with the case r = 3/5 = 0.6 and our 
fitting attempts are summarized in figure HI Figure IU(a) shows a simultaneous fitting - 
including all data points - using the shift behavior (JSj) for the estimation of the critical 
temperature and the correlation length exponent. Now, we have used the pseudocritical 
temperatures T\z\* av of the specific heat (Z = C) and magnetic susceptibility (Z — 
As shown in the panel, the estimated value for the critical temperature is T c = 2.2245(7) 
in excellent agreement with the exact value. Respectively, the estimation of the inverse 
of the correlation length exponent is \ jv = 1.01(1), which provides a value for v of the 
order v = 0.99(1), very close to the value u — 1 of the pure model. 

In panel (b) of figure H] we illustrate the data of the averaged specific heat at the 
exact critical temperature. The solid and dotted lines are corresponding fits of the form 

[CV(T = T c ) ~ C*! + C 2 In (In L) (9) 

and 

{C] av (T = T c )~C oc + bL a /». (10) 

As it clear from this graph, one may not easily discern between the two lines, although 
a more careful analysis indicates that the double logarithmic function is more stable. 
We have performed both kinds of fits for three sets of data points L min — L max . We 
fixed L max = 120 and varied L min as follows: L min = 20, L min = 50, and L min = 80. 
We have observed that although the values of x 2 /DoF are comparable between the 
two functions ([9]) and (fTUI) . the coefficient C 2 of equation ([9]) seems to be stable 
(C*2 — 1.67(3)), whereas the exponent a/u of equation (TTDT) fluctuates, with increasing 
L min , in the range -0.3 to -0.2: f (L min = 20) = -0.24(2), %L min = 50) = -0.21(4), 
and ^(L min = 80) = —0.30(3). Although for this case we can not conclusively 
discriminate between the two alternatives, the inset in panel (b) shows that the specific 
heat data fits well to the double logarithm for lattice sizes L > 60. Finally, in panels 
(c) and (d) we plot the average magnetization and the average magnetic susceptibility 
at the critical temperature respectively, as a function of the lattice size L in a log-log 
scale. In both cases a linear fit is applied for L > 20 giving the values (3/v — 0.125(3) 
and 7/V = 1.754(6), i.e. the ratios of the magnetic exponents (3/v and 7/V maintain 
the values /3 p /z/ p = 0.125 and 7 P /V P = 1.75 of the pure model. 

We proceed to present our results for the strong disorder strength r = 1/7 = 0.142. 
As was already stated above, for this value of r we have extended our simulations up 
to sizes L = 200. Figure [5] summarizes now the critical behavior for the case r = 0.142 
and, in particular, in figure [5]^a) we estimate by a simultaneous fitting attempt for the 
larger lattices (L > 60) the critical temperature and the correlation length exponent. 
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Figure 5. The same with figure [4] for r = 1/7 = 0.142. Now data for lattice sizes up 
to L = 200 are presented. All fitting attempts shown are performed for L > 60. 



The estimated value for the critical temperature is T c = 1.7816(30) in good agreement 
with the exact value. The production of the above accurate estimates for the critical 
temperatures, even in this strong disorder case, constitutes a concrete reliability test, 
in favor of the accuracy of our numerical scheme. The estimation of the inverse of 
the correlation length exponent is \jv — 0.98(3), providing a value for v of the order 
v = 1.02(3), which within error bars agrees with the pure model's correlation length 
exponent value. In panel (b) of figure [5] we plot the data of the averaged specific heat at 
the exact critical temperature. Again, the solid and dotted lines are corresponding fits 
of the forms ([9]) and ( flOi) for the larger sizes (L > 60). Our analysis for the quality of the 
fits for the three sets of data points L min — L max , with L max = 200 and L min = 20, 50, 
and 80 indicated a good trend for the values of x 2 /DoF for the double logarithmic fits ([9]) 
in the range: 0.1 — 0.3. Further reliability in favor of the logarithmic corrections scenario 
is provided by the stability of the coefficient C<i- C2{L min = 20) = 0.412(4), C2(L min = 
50) = 0.410(6), C 2 (L min = 80) = 0.418(7), whereas the value of the exponent a/u of 
the power law (fTOl) approaches zero with increasing L min : ^(L min = 20) = —0.27(4), 
%L min = 50) = -0.20(6), and f (L min = 80) = -0.02(3). From the above, we may 
conclude that our numerical data are more properly described, at least for the larger 
lattice sizes studied, by the double logarithmic form ([9]). The inset in panel (b) of 
figure [5] shows again the specific heat data as a function of the double logarithm of 
the lattice size for L > 60. The solid line shown is an excellent linear fit. It should 
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be noted that, the above power law behavior (a/u — > 0_) is in full agreement with 
the earlier observation of Wang et al [25] in the strong disorder regime (r = 1/4 and 
r = 1/10) [25]). Note also that, analogous results have been presented in figure 2 of 
reference [39] for the site diluted Ising model, where sizes up to L = 256 have been 
considered. Finally, in panels (c) and (d) we plot the average magnetization and the 
average magnetic susceptibility at the critical temperature respectively, as a function of 
the lattice size L in a log-log scale. In both cases, the solid lines shown are linear fits 
for L > 60 giving again within error bars the values of the pure model, j3/u — 0.125(6) 
and j/u = 1.752(9), respectively. 

Summarizing, we may point out that the difficulties observed in the case of weak 
disorder in discriminating between the two scenarios have been surpassed by considering 
the strong disorder case and larger lattice sizes. In particular, we take as evidence in 
favor of the picture emerged from the theoretical work of Dotsenko and Dotsenko [H] 
and the improved versions of Shalaev [19], Shankar [2T] . and Ludwig [22], the stability of 
our fitting attempts using the double logarithmic law. Similar conclusions have also been 
reported by the extensive Monte Carlo studies of Wang et al [25] for the random bond 
model and also of Selke et al |43j . Ballesteros et al [39], and Tomita and Okabe [IB] for 
the site diluted model, by the transfer matrix approach of Arao Reis et al [361 HO] for the 
random bond case, and very recently by Kenna and Ruiz-Lorenzo [58] via an alternative 
approach that involves the density of Lee- Yang zeros of the site diluted model. 

5. Conclusions 

The effects induced by the presence of quenched bond randomness on the critical 
behavior of two Ising spin models in 2D have been investigated by a sophisticated 
entropic scheme based on the Wang-Landau algorithm. For the random bond square 
SAF model we have extracted accurate estimates for all critical exponents and two 
values of the disorder strength. These values verify hyperscaling, satisfy the Chayes 
et al inequality [5], and obey very well the weak universality scenario for disordered 
systems |113j . Furthermore, the strong saturating behavior of the specific heat clearly 
distinguishes this case of competing interactions from other 2D random bond systems 
studied previously. For the marginal case of the random bond Ising model, our findings 
favor the well-known double logarithmic scaling scenario and suggest that the pure 
system behavior, v — 1, is recovered in the asymptotic limit. Here, the estimated 
critical temperatures, in both cases of disorder, are in excellent agreement with the 
exact values obtained by duality reflecting the accuracy of the implemented entropic 
scheme. Encouraged by this latter observation, we are currently carrying out a similar 
study of bond disorder effects in 2D systems undergoing first-order phase transitions. 
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